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solving scalar conservation laws. By making suitable approximations we obtain 
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1 . Introduction 


In solving hyperbolic conservation laws of the form 

(1.1) + f(u)^ = 0, 

it is desirable to have a method that is at least second-order accurate in 
smooth regions of the flow and that also gives sharp resolution of 
discontinuities with no spurious oscillations. 

There are two basic approaches which have been used to derive difference 
schemes for these problems. The first is purely algebraic. A scheme is 
defined by its coefficients or numerical fluxes, and algebraic relations are 
derived which guarantee that certain desirable properties hold. 

The second approach is more geometrical, in that the structure of certain 
special solutions to (1.1) is heavily used. The classic method of this type 
is Godunov's method [3], based on the solution to Riemann problems. Van Leer 
[10], [11] with his MUSCL schemes, generalized this method to second-order 
accuracy by using discontinuous piecewise linear approximations at each time 
step. Higher order geometrical methods, such as the piecewise parabolic 
method [1], have also been used. 

Recently, however, most theoretical progress toward deriving second-order 
schemes with desirable properties has been made using the algebraic approach. 
Harten [4] introduced the concept of a total variation diminishing (TVD) 
scheme which guarantees that monotone profiles remain monotone. He derived 
conditions on the coefficients of a scheme that guarantee the TVD property and 
second-order accuracy. Since then many other second-order accurate TVD 
schemes have been constructed, e.g., [13], [15], [17], [18], but always using 


algebraic methods. 
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Our goal here is to derive a simple second-order accurate TVD scheme 
using a geometric approach in the spirit of Van Leer's work. This is 
accomplished by choosing the slopes properly in the piecewise linear 
approximation and then by also approximating the flux function f in (1.1) by 
a piecewise linear function. The purpose of the latter approximation is that 
the modified equation with piecewise linear initial data can be efficiently 
solved analytically. 

The method we derive is written in conservation form and can be viewed, 
algebraically, as a "limited" version of the Lax-Wendroff method. For a 
linear problem it agrees with one of the flux-limiter methods studied by Roe 
[17] and Sweby [18] but differs for nonlinear problems. 

One advantage, we feel, of the geometric approach is that it gives more 
insight into the behavior of algorithms. It may make it easier to show, for 
example, that the resulting numerical solution satisfies the entropy 
condition. Toward this end we choose a geometric form of the entropy 
condition, namely that solutions satisfy the spreading estimate 

(1.2) - u(y,t) £ ^ If x>y 

X y t 

for some constant c. Oleinik [12] has shown that weak solutions to (1.1) 
satisfying (1.2) are unique. 

In Section 3 we prove (1.2) for the approximations produced by Godunov's 
method. The analysis shows, moreover, that these grid functions satisfy 
(1.2) with the correct constant c away from sonic points and points where 
the CFL condition is binding. Thus, Godunov's method spreads rarefaction 
waves at the physically correct rate most of the time. At the sonic point. 



This 


(1.2) is satisfied with a constant that is two to four times larger, 
causes a sonic rarefaction to develop a kink, or "dog leg" at the sonic point, 
as has frequently been observed without explanation in computations. 

Unfortunately, the technical arguments used to prove (1.2) for Godunov s 
method do not carry over directly to the second-order scheme so that we have 
not obtained the spreading estimate in this case. However, by considering the 
sonic rarefaction case we will argue that entropy-violating shocks cannot 
persist. Moreover, numerical results look very good with spreading at the 
correct rate everywhere, including at the sonic point. 


2. Godunov's Method and Second-Order Extensions 

We consider the scalar version of equation (1.1) and will always assume 

that the flux function f is convex: f" > 0. We will denote the numerical 

approximation to the solution u(xj,tj^) by Uj. Here Xj = jh and 

t = nk where h and k are the mesh width and time step, respectively, 
n 

Since we will be discussing formulas for a single step from t^^ to we 

will generally drop the superscripts and replace U? and by Uj and 

Uj , respectively. 

To take a single step with Godunov's method, a piecewise constant 
function w(x,t^^) is defined which takes the value Uj in the interval 
Ij = (Xj_ I /2 equation (1.1) is then solved exactly up to time 

tn+i with this initial data to obtain w(x,tjj+p. This can be done easily 

if k is sufficiently small by solving a sequence of Riemann problems. The 
new approximation IJj is then obtained by averaging the solution w(x,tjj+p 


over Ij s 
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( 2 . 1 ) 


Vi/ 


V2 


w(x,tn^l)dx. 


X. 1, 

J- V2 


The method can also be rewritten in conservation form as 


( 2 . 2 ) 


= Uj - I [F(U;j) - F(U;j-l)J 


where 


(2.3) 


^n+1 


F(U 


/ UTl 

f(w(x^^l, ,t))dt 


n 


is the "numerical flux" across Xj_^ , and depends only and . For 
a scalar conservation law this expression for the flux can be simplified. 
Assuming for convenience that we are away from the sonic point (so f'(u) ^ 0 


(2.4) 


F(U;i) = 


) we have 



f(u^) 

if 

f' > 0, 


if 

f' < 0. 


The formula for the sonic case is only slightly more complicated. Using this 
formulation allows Godunov's method to be applied with any size time step for 
which the Courant number is less than 1 . For a scalar conservation law the 
Courant number v is defined by 


V = max |f'(u)[. 

For V < 1, Godunov's method is TVD, i.e., 
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TV(U) < TV<U) 


where the total variation is defined by 

TV(U) =2 lUj+i - Ujl. 

j 

In modifying Godunov's method to obtain second-order accuracy, we follow 
Van Leer [10], [11] and replace the piecewise constant function w(x,tj^) 

defined above by a piecewise linear function which we will denote by v. This 
replacement is conservative provided v is of the form 


(2.6) v(x,t ) = U. + s.(x - X ) for x £ I . . 

J J J J 

We wish to pick the slopes Sj so that the total variation of v is the same 

as with s. =0. One simple choice is 
J 

/o if -UjKDj -Uj.j) <0 

(2.7) s. = < 

^ / U . - u, u. - u._. 

VsgnCUj^^ - Uj) min{lJ — ^ |-J — } otherwise. 

Using this choice of slopes we can define the following algorithm for solving 

( 1 . 1 ): 


Algorithm 2.1 . 

1) Determine v(x,t^) based on {Uj} using (2.7). 

2) Solve (1.1) exactly with initial data v(x,t^j) to obtain v(x,tj^+].^* 

3) Average v(x,tjj^j) as in (2.1) to obtain Uj . 
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Each of these steps is total variation diminishing; so this defines a method 
which is TVD. Moreover, it can be shown that this method is second-order 
accurate in smooth regions, at least away from extreme points of u. 

Unfortunately, this is not a practical method in most situations since it 
requires solving the conservation law (1,1) exactly with piecewise linear 
initial data. This is more difficult than solving Riemann problems. Various 

modifications can be made to Algorithm 2.1 to give a more readily Implemented 
method. 

Here we introduce a variant which remains second-order accurate and TVD 
and is easily Implemented. We solve (1.1) with the piecewise linear initial 
data (2.6), but only after modifying the flux f in (1.1) to make this 
tractable. Specifically, we replace f by a piecewise linear function in 
such a fashion that computing the flux across Xj+ reduces to solving a 
linear problem with piecewise linear data. The solution to this problem is 
easy to derive. 

To compute the flux across , we first compute the slopes (2.7) and 
consider the function v(x,t^j) In (2.6). Since the sonic point causes 
difficulties, we delay discussion of this case to the end of this section and 
begin by assuming that f'(v(x,t^)) ^0 for x € Ij U . Recall that we 
are always assuming f is convex. 

Set 

^i “ ^i * 7 ^®i* 


By virtue of our choice of slopes (2.7), the points uj, Uj, uJ^j, are 

monotonically ordered (though two or more may coincide). Let g(u) be a 


piecewise linear function which interpolates f(u) at these four points and 


set 


(2.9) 


f(U^) - f(U^) 

u^-u- 




if s^ * 0 


if Sj^ = 0 


for i = j, j+1 , so g' is the slope of g(u) between and U^. 

Now consider the problem 


( 2 . 10 ) 


Vj. + g(v)jj = 0 


for t < t < t with the piecewise linear initial data v(x,t_). We can 
easily compute the flux across during this problem, 
which we will denote G(U;j)* Since we can rewrite (2.10) as + g^(v)v^ = 0 
and constant for £ t ^n+1’ that 


( 2 . 11 ) 


- <t - t_)sj g: 

v(*j+V2.t) - j 

(“j+l ■ ®j+l 


if f' > 0 


if f' < 0. 


Furthermore , 


(2.12) g(u) = 


fCuJ) . (u - u;>gj 


£(Vi> - <» - Vi)8j'.l 


for u £ lnt[U, ,ut] 




for u 
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where int[a,b] denotes the interval with endpoints 
and these formulas, we' find that the flux is 

a and 

b. Using 


in+1 

G(U;j)=iy g(v(x^^y^ ,t)]dt 

t 

n 



(2.13a) 


if f' 

> 0, 

(2.13b) 


if f' 

< 0. 


Notice that if = 0 for all j, i.e., if we use piecewise constant initial 
we recover the flux (2»4) of Godunov s method* Also note the similarity 
of both (2.13a) and (2.13b) to the flux of the Lax-Wendroff method which can 
be written in the form (2.2) with flux 

In fact, for smooth solutions, 

(2.15) G(U;j) = FLy(U;j) + O(h^) 

2 

and the 0(h ) term is a smooth function of j (except where Sj = 0, i.e., 
at extreme points of U). It follows that in computing 

(2.16) [^(U;J) - G(U;j-D], 

the 0(h2) terms cancel to 0(h2) showing that our scheme agrees with Lax- 


(2.14) =i(f(U.) + f(U._^,)} -i 


3+1 ' 
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Wendroff to O(h^) locally and hence is second-order accurate in smooth 
regions, except near extrema of U (where it seems that all known second- 
order TVD schemes reduce to first order accuracy [15]). 

The use of the limited values Sj and gj in (2.13) rather than the 
corresponding expressions in (2.14) gives us a method that, unlike Lax- 
Wendroff, is TVD. To see that this is so, note that it suffices to check the 
following conditions: 


(A) If Uj is a local maximum, i.e., ^j-1 — ^j+1 

local minimum), then ^ (resp. 2. 

(B) If 1 1 ^j+l ^j-1 ^ ^ j ^ ^j+1^ then 

Vi < Uj < Vi > "j > "1+1>- 

We are still assuming that we are away from the sonic point, specifically that 
f^(u) has one sign on I. The sonic case is discussed below. 

So suppose, for example, that f"" > 0 and that ^j-1 other cases 

are completely analogous). 


If Uj > U.., then we must check (A). In this case s^ = 0 and 
j — J+1 


G(U 


;j) = f(Uj) while Sj-i ^0 


and 


G(U;j-l) = f(Uj_p - Ik 1 


So, 


Uj = Uj - ^ [G(U;j) - G(U;j-l) j 

< U. 

- J 


and (A) is satisfied 



If Uj < then we must check (B), We require, of course, that the 

Courant number be less than 1; hence Uj depends only on values of v(x,tjj) 
In Ij_^ Ij ^j+1* Since Uj_j £ U j £ we have 0 for i = j-1, 

j, j+1; and 


(2.17) 


Vi < “j-i < «j < »; < Vi 1 


The new value is determined by averaging the exact solution to (2,10). 
In our derivation we defined the piecewise linear flux g(u) locally; it had 
one definition in computing G(U;j) and a different definition in computing 
G(U;j-l). However, by the condition (2.17) these definitions are consistent 
in the region where they overlap, and so we can define a single function 
g(u) to compute both fluxes. Specifically, we can take 


+ (u - 




g(u) = <( f(Uj) + (u ~ 

f(u!^l) + (u - U+^^)g'^^ 


for u < U. , 

- J-1 

for U* < u < U* 

J-1 - - J 

for U. < u, 

J - 


where 


u* . (f(u*) - £(u;^j) + - g;u;)/(gj^, - gp 


for i = j-1, j are the points of intersection of the piecewise linear 
By the convexity of f and (2.17), we find that 


segments. 
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so that this definition is consistent with the previous definitions and gives 
the correct fluxes. Moreover, since < - ®j+l ’ this piecewise linear 

approximation is also convex. It follows that is obtained from 

i s j-ij j, j+1, by solving a conservation law with a convex flux function 
and, hence, by the TVD property of the exact solution. 




so condition (B) is satisfied and the method is TVD. 

We now turn to the sonic case and derive formulas for the flux across 
Xj_j_ jy when the sonic point uq lies in int[Uj j . In some cases the 

previous formulas (2.13) are still valid and, to avoid repeating these 
expressions, we define 


(2.18) 


G, = f(U"^) - i ks (g' )^, 

J J ^ J J 

®j+l "" " T ^®j+l^®j+l^ * 


We first note that if g' and g'^^ have the same sign then the 
linearized problem (2.10) can be solved just as before and the flux agrees 
with (2.13), 



if g' > 0, g'^^ > 0 
if gj < 0, g'^^ < 0. 


G(U;j) 


- 12 - 


If g' and have different signs, then we must be more careful. 

If gj < 0 and g'^j > 0, then the discontinuity at Xj^. is a sonic 

rarefaction. By solving (2.10) in this case we find that the flux is f(ut) 
if Uj Uq £ Uj and f(Uj_j_j) if £ '^0 — ^j+1* remaining case, 

^ ^0 ^ ^j+ 1 * sonic point lies within the discontinuity in v(x,t ). 

In this case we must take special precautions to ensure that the rarefaction 
wave spreads properly and an entropy violating shock does not persist (see 
Section 3). Rather than using the usual piecewise linear flux g(u) we 
include another Interpolating point (uq, f(uQ)) in g(u). The flux is then 
simply f(uQ). 

These last three expressions for the flux can be conveniently combined to 

give 

G(U;j) = f(v^) if g' < 0, g'^j > 0 


where 


(2.19) 


Vq = min[max(Uj,UQ), 


Now suppose gJ < 0 and g'^^ >0. In order to solve the linearized 
problem (2.10) we must also specify g(u) for <. u ^ in this case. 
We take another linear segment interpolating f at these points with slope 


(2.20) 


"j+ Vz 


(f(Uj^j) - f(U+)' 

"j+1 “ 

f'(uj) 


“ Vi ‘ ”r 


« Vi-V 
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In the solution to (2.10) the discontinuity at \j^ is now a shock which 

may propagate to the left or right or be stationary. If the shock moves to 
the left (for all t C (tjj.t^^^l)) , then the flux across 1/^ is Gj+i 

while if it moves to the right the flux is G j . Unfortunately, since the 

initial data is not constant on the two sides of the discontinuity, the shock 
may switch direction and cross Xj+ at some time t C which 

point the flux is discontinuous. This is most easily visualized by first 

considering the multi-valued solution obtained in solving the linearized 
problem- and then inserting a shock according to the equal area rule. By 
taking the slopes Sj and Sj^j quite different one can construct examples 
where the shock is first on one side of Xj^ 1^2 and later on the other. 

For simplicity we ignore this possibility and always use Gj or Gj+j 
depending on the initial motion of the shock. This is the one situation in 
which we do not use the exact solution to the linearized problem, but 
experimentally this approximation seems to work well. 

If v(x,t^) is discontinuous at *j+ 1/2 » ^ ^j+1’ ^1*®” 

initial motion of the shock is determined by the motion of the discontinuity 
in the multi-valued solution and hence by the sign of g'^ . We use 


o(Oij) - ■! 

( Vl 


if 


If g' 1, =0 then the discontinuity is stationary and the initial motion of 

j+ V2 

2 ^ 1 

the shock is determined by the relative sizes of aad Sj^^(gj^j) . 
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Specifically, 

If SjCgp' > Sj^iCgp,)^ 
otherwise. 

It turns out that this is also the appropriate formula if ut = • In 

this case sj = Sj+j and we are simply comparing the magnitudes of g: and 

®j+l • 

We have now derived formulas for every possible case. Luckily all of 
these formulas can be summarized quite neatly as follows: 

Algorithm 2.2 

Uj = Uj [G(U;j) - G(U;j-l)J 

where 


0(U;J) = 


J+1 


» If gj>0, SPV 2 <V, and 


gPj < 0: 


h 

G(U;j) = / 

(®j+l 

2) Otherwise: 

T" 

G(U,J) . n 

(f(Vo) 


If apgp2 > 

Otherwise. 


if 

gj > 0. 

®j+ V 2 

if 

gJ+1/2 1 

®j+l - 

if 

gj < 0. 

> 
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Here we have used the expressions (2.7), (2*8)) (2»9), (2«18), (2«19), 
and (2.20). These formulas cover all cases; the sonic shock, sonic 
rarefaction and also the usual nonsonic case (2.13). However, in implementing 


this method it is of course best to use (2.13) whenever g' g'^^ >0. We 

should compute g'^ y and perform the various tests above only in the 
relatively rare sonic case. 

It is possible to show that the method remains TVD even near sonic points 
when these formulas are used. This is done in precisely the same way as 
before but is slightly more complicated since several cases must be 


considered. We omit the details. 

Numerical experiments confirm that the method is second— order accurate 
and TVD. To check the second-order accuracy we applied the method to Burgers' 
equation u^ + uUj^ = 0 with smooth (sine wave) initial data and periodic 
boundary conditions. Both the Lj and norm of the errors decrease at 

the correct rate as the mesh is refined. 

Figures 2.1a and 2.1b show the results of a typical calculation in which 
a shock forms. Again the method is applied to Burgers' equation with initial 


data 



X < .5 
X > .5 


and periodic boundary conditions. The discontinuity at x = .5 spreads into 
a rarefaction fan, and the smooth decreasing profile sharpens into a shock. 
Figures 2.1a and 2.1b show the results at time t = 0.2 and t = 0.4, 
respectively. For comparison. Figure 2.2 shows the results of Godunov's 
method. Notice the improved accuracy in the smooth portion of the solution 
with the second— order method and the lack of oscillations near the shock. 
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Godunov's method suffers in particular from a lack of smoothness in the 
rarefaction wave at the sonic point Ug = 0 which does not occur with the 
second-order scheme. This is discussed in Section 3. 

It is interesting to compare this method to the flux— limiter methods. We 
find that for a linear problem it is the same as one of the flux-limiter 
methods of [18] but that it differs for nonlinear problems. 

First consider the linear problem 
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where 


(2.23) 


^(r) = 



Then (2.21) becomes 


if r<0, 
if 0 < r < 1 , 
if r > 1. 


(2.24) “ '’^^j+1 " ■ T 


where A o) = u - u. ,. This is precisely the flux-limiter method of [18] 
- j 3 J-1 

with the so-called "minmod" limiter given by (2.23). We note in passing that 
other flux-limiter methods are given by different choices of the limiter ^ 
and that taking (j)(r) = 1 for all r gives the Lax-Wendroff method. 

The flux-limiter method is extended to the nonlinear problem (1.1) by 
generalizing (2.24) to 


(2.25) “i “ “j ■ ¥ 




in the case f' > 0, where 


and 


'’ 3 + >/2 


■«Vl> • 

Vl-"3 . 


<1 - 

■^3 “ a - - f(Uj)) • 


By contrast, our method in the same situation gives 
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"j =«J -iA.[f(uJ) 4„Sj(gp2j. 

Using the definitions of the various quantities appearing here, we can 
rearrange this to obtain a form similar to (2.25): 

(2.26) T. = Uj - ^ (£(0j) - £(up) - ‘ ^ 4jU - I 8p(£(«J) - £(up)J. 

This is very similar to (2.25) but the limiting is done in a different manner. 


III. Spreading of Rarefaction Waves 

Weak solutions to conservation laws are not necessarily unique. In 
general there is some additional condition, such an an entropy condition, 
required to identify the unique physically relevant solution [5], [9]. For 
the scalar conservation law (1.1) with a convex flux f, such conditions are 
well known in several equivalent forms. One form considered by Oleinik [12] 
requires that the solution satisfy the spreading estimate 

(3*1) u(x,t) - u(y,t) < — y . 

— at 

for all X > y and t > 0 where a > 0 is some constant. In fact, one can 
take a = a where 

(3*2) a = inf f”(u) 


and a > 0 by convexity. We can define a locally to obtain more precise 
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information on the rate of spreading of rarefaction waves in different regions 
of the solution. Note that at points where u,j(x,t) exists we obtain 

(3.3) u (x,t) < , 

We would like to prove an estimate analogous to (3.1) for the numerical 
solutions generated by a particular scheme as the mesh is refined with k/h 
held fixed. If there exists a constant c > 0 such that 


(3.4) 


- u? < 

j+1 J — cn 


for each point on every grid, then the limit solution satisfies (3.1) with 
a = ch/k and hence satisfies the entropy condition. The possibility of 
proving estimates of this form was independently noticed by Tadmor [19]. He 
proved the estimate (3.4) for the Lax-Friedrichs scheme using essentially the 
same technique. 

This form of the entropy condition seems easiest to deal with when 
studying second-order schemes of the type considered here. Moreover, by 
obtaining an estimate of the form (3.4) we can compare the rate of spreading 
in the numerical solution with the correct rate. Ideally we would like 
c = ak/h in (3.4) . 

Our interest in obtaining such quantitative information stems from the 
observation that rarefaction waves computed with some numerical methods, 
including Godunov's method, do not always spread at the proper rate in spite 
of the fact that the entropy condition is satisfied. This difficulty is most 
frequently observed at the sonic point. As we will show with Godunov's 
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method, rarefaction waves spread with at best one-half the correct rate in 
this region. This leads to a kink in the rarefaction wave at the sonic 

point. This is frequently observed in practice and has been termed a "dog- 
leg" by Sweby [18]. For an example see Figure 2.2 where we have applied 
Godunov's method to the Burgers equation u^ + uu^^ = 0. 

Although our main interest is in the second-order methods of Section 2, 
we will begin by analyzing Godunov's method in some detail. This will provide 
a basis of comparison and also provides some insight and quantitative 
information on the dog— leg phenomenon. 

Let 

d" = (u" , - u’?). 

J ^ J+1 

Then our goal is to find a constant c > 0 so that 

(3.5) d" < — 

— cn 

for all n. In fact we will consider only a single time step and, as in 
Section 2, replace d" and by D and ^ respectively. We will 

determine a constant c > 0 for which 

(3.6) D _< D - cD^ 

from which (3.5) follows by induction. 

We also let 



U. 


(3.7) 
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and note that 


(3.8) 


D. 

J — ko 


for all time steps by virtue of the Courant number restriction 


1^ f '(u) I < 1. Here a is as in (3.2) and (3.8) follows from 

U. 

k . k f 

h”3 


U.J.1 

• 1+1 


f"(Od5 




< 2 . 


Now consider the mesh point and suppose to begin with that 

f'(u) >0 on LJ Ij Then we will show that 


(3.9) 


D. < D - cD 
J - 


with c ak/2h and that, in fact, we can generally use c = odc/h. Since 
f' > 0, applying Godunov's method (2.2) and (2.4), 


(3.10) 


■ Vi "h 




Subtracting these and using (3.7) gives 
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" “j - H ^ I ”j ^ J “j-l 


1 ^2 


where we have expanded both f(Uj^j) and fCUj^j) in Taylor series about 
Uj and Uj_j^ £ ^ ^ <. Rearranging gives 


(3.11) 


“ 2 h ['*j '*' °j-l 


Since ^ ^ f'(Uj) ^ 1, the first two terms are a convex combination of Dj 

and Dj_i and consequently bounded by 


(3.12) = max(Dj_j^ ,Dj ) . 

Dropping the term corresponding to the smaller of Dj and from the sum 

in brackets in (3.11) and using (3.2) gives 

(3.13) Dj < D* -i^ aD^ . 

The term on the right of this inequality is an increasing function of D* 

for D* satisfying (3.8), and so we can replace D* by D = max Dj and 

conclude that (3.9) holds with c = ak/2h. 

Moreover, except near the extreme points of U (l.e., the edge of the 
raref unction wave) we can do better than this. Typically in the interior of 
the rarefaction wave we have • 
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If in fact 

(3.14) = Dj , 

then by retaining both terms in the sum in brackets in (3.11) we remove the 
factor V 2 ill (3.13) and obtain 

(3.15) D. < D - ^ a D^. 

Even if (3.14) does not hold exactly we can consider modifying the data by 
increasing or decreasing so that (3.14) does hold. Then (3.15) 

holds for the modified D ^ , but it is easy to verify (using the fact that 
Godunov's method is monotone) that the original is bounded above by the 

modified so (3.15) holds for the original data. This argument works 

provided modifying the data does not violate the Courant number restriction, 
which it might near extreme points of U. 

Exactly the same arguments can be applied to the case where f'(u) < 0 
on IjU conclude that, away from the sonic point, 

rarefaction waves spread at the correct rate except perhaps near the edge of 
the rarefaction wave if the Courant number restriction is binding there. 

Now consider the sonic case, where j ^ ^ ^'+1* 

Godunov's method in this case gives 

Vi ■ Vi - K 


Applying 



-24- 


Subtracting gives 

(3.16) Dj = Dj - I F 2 

where 

Since f'CuQ) = 0, expanding f(Uj^j) and f(Uj) about uq gives 

(3-17) ^2- I 1*^1 

where A... = U.., - u„, A. = u„ - U., and U. < C. < u» < 5..i < U.,, . 

j+1 j+1 O’j 0 3* J — J— 0 — ^3+1 — 3+1 

Since A. + A.,, = D. we have 
J J+1 J 

(3,18) -i- a Dj < F 2 < j a Dj. 

So from (3.16) we obtain (3.6) with a value of c which is at worst V4 the 
correct value and at best one-half the correct value. The worst case occurs 
when Aj = ^j+1 ” T ^ j * ^*®*» 1-^ *-1^® sonic point falls half way between Uj 
and Uj+i* Th® best case occurs when A^ ® either Uj or Uj+i 

is equal to uq. 

Summarizing these results, we can obtain a global bound of the form (3.5) 
with c = ok/4h which shows that the entropy condition is satisfied. 
Moreover, we generally have spreading at the correct rate except near the 
sonic point, where the rate is at best one-half the correct rate. 
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These results can be seen more geometrically by considering Figures 3.1 - 
3.4, which for simplicity we have shown for Burgers' equation where 
rarefaction waves are linear in x. Moreover, we have taken ^ j • 

Figure 3.1 shows the initial conditions U as the solid piecewise constant 
function. After solving Burgers' equation with this data, we obtain the 
dashed line. It is this function which is averaged to give , indicated by 
dots. Clearly Uj and both decrease, but decreases by a greater 

amount so that the difference Dj also decreases. A little reflection shows 
that Dj decreases by 1/h times the area of the shaded rectangle drawn in 
Figure 3.2. The height of this rectangle is D j , and the length is easily 

computed to be kD^a, where a = f" = 1 for Burgers' equation. So, 


D. = D. - ^ a (D.)^ 

J J h 2 


which agrees with (3.15). 

Now consider the sonic case and suppose first that Uj = 0. Then we have 
the situation in Figure 3.3, and Dj decreases by 1/h times the shaded 
triangle, which is one-half the area of the rectangle in Figure 3.2. This 
accounts for the spreading rate being one-half the correct rate in this case. 

Finally, considering the worst possible sonic case, where 
Uq = -j (Uj + we obtain Figure 3.4. Again, Dj decreases by 1/h times 

the shaded area, which is now one quarter the original area. 

We now turn to the second-order accurate scheme given by Algorithm 2.2. 
Unfortunately, we have not yet been able to prove the desired spreading 
estimates in general. The approach used above for Godunov's method does not 
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apply directly since it is possible to construct initial data for which the 
difference actually increases in a single time step. This undesirable 
behavior does not persist in later time steps, and experimental evidence 
indicates that rarefaction waves do spread at the correct rate asymptotically, 
but clearly we can no longer use a bound of the form (3.9) to prove this. 

Nonetheless, we can gain some theoretical confidence in the method by 
considering the sonic rarefaction case. In practice this is the most 
worrisome case; violation of the entropy condition is usually manifested in 
the form of a sonic shock which fails to spread. Intuitively, we expect that 
the second-order method, being an extension of Godunov's method, will not 
permit such behavior. In fact, if we consider data for which 


so that the sonic point lies within the discontinuity between Uj and 
we find that 

- £(u]^j) - i k 

(3.19) G(U;j) = f(uQ), 

G(U;j-l) = f(uj) - i k Sj(gJ)^ > f(Uj). 


The inequalities here follow from the Courant number restriction and the 
convexity of f. For example, 0 > kg' ^ -h and by the definition of g'. 


f(Uj) = f(u") + hSj g' 
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so that 


f(Uj) < y (£<Uj) + fWj)) 


- *2 Sj' 

< £(UJ) -|kSj(g')2 


and similarly for the other inequality. 
From (3.19) we obtain 


= ^j+1 " ■ 2G(U;j) + G(U;j-l)J 

> Vi - "j - K («Vi> ■ 

Comparing this to (3.16) shows that the second-order method has spreading with 
at least the same rate as Godunov's method in this situation. 

In practice, it seems to be much better than Godunov's method as in 
Figure 2.1. The prominent kinks in the rarefaction wave computed with 
Godunov's method are entirely missing in the second-order calculation. 

As a final comment about entropy, we note that by placing additional 
constraints on the slopes Sj we could obtain the standard entropy inequality 
for all entropies of the form ri(u) *= |u - c| with c any real constant. 
Kruzkov [7] has shown that this is sufficient to guarantee uniqueness. We 

will sketch this only briefly since we do not feel that such a modification is 
necessary in practice. 

For Godunov's method it follows from Jensen's inequality for convex 
functions that the entropy condition is satisfied. If (h»q) is any convex 
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entropy pair (see, e.g., [9]) then 

J h(w(x,t^^p)dx 

1/2 " j - 1/2 

where, as in Section 1, w(x,t) is the exact solution to (1.1) on (^n»*'n+l^ 
with piecewise constant initial data Uj . It follows that 


(3.20) n(U^) - n(Uj) < ^ 




which is the discrete form of the entropy inequality. 

For the second-order method we must also consider the step going from 
to the new piecewise linear function v(x,t^^P ~ lij ^ 

(Xj_ , Xj^y^). For the Kruzkov entropies we would like to show that 


(3.21) 


P+V2 

J I"’'- 

"j- V2 


"n+l’ - <'1'“’' 




1/2 


|v(x,tn^l) - c|dx 


"i- 1/2 


since we can then show the discrete entropy inequality 


Ej - Ej < ^ [Q(U;j+D - Q(U;j)J 


where 
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^j- 1/2 


Q(U;j) 



,t))dt. 


Here v(x,t) is the solution to the linearized form of (1.1) with piecewise 
linear initial data. 

However, (3.21) may fail to hold for some values of c if the slopes 
Sj defining v are chosen according to (2.7). Since (3.21) does hold if 
Sj = 0, it should also hold for Sj sufficiently small. In particular, we 
can show that a sufficient condition for (3.21) is 


(3.22) 


I=ii <*1V 


Since in smooth regions of the flow we expect s^ » this is not much of a 
restriction and the modified method should still be second-order accurate. On 
the other hand, the restriction (3.22) would be difficult to impose in 
practice and does not seem necessary since we have not encountered any 


difficulties with the unmodified method 
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Figure 2.2 Godunov's method 
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